C.................... RSMNSD.FOR ..................................... 
C      PROGRAM FLIP
C.... this main module for the Field Line Interhemispheric Plasma (FLIP)
C.... model. It reads initial conditions and sets up the
C.... coordinate system and neutral atmosphere, note s(4500) for d300
C
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      INTEGER ION,NEQ,WUNIT8,IRIHED,ICAUTN,IVERT,IVPERP,NCOL_EXB,JDAY
     > ,UTORLT,NOCON,AUR_FILE,APFILE,JZLO,JJ,FSTART,JZ,JBLATPR,J250,I
     > ,FEXISTS
      INTEGER NUMUNIT,UNITNUMS(22)     !.. Unit numbers for writing rate warning
      REAL SEC,DEC,ETRAN,F107,F107A,AP,BLON,D,T,CHI,SATN,SATF,SAT,EUV
     > ,UVFAC,GLON,GLONN,GLONF,CHIN,CHIF,GLATN,GLATF,GLATJ,AP3HR(33001)
     > ,GRD,RATFAC,SW(25),UTHRS,UHED(2),SATEQ,L_STAG,L_SHORT
     > ,RE,PLAT,PLON,PI,ANGVEL,TIMEON,TIMEOF,AURALT,ZTOPTE,L_PHASE,SX
     > ,TD_PHASE,DELTIM,PHASE,SECWF8,GLOND,GLATD,BCOMPU,KPREAL,ZFLUX(37)
     > ,SATX,CHIX,GLONX,GLATX,SATNSAV
      !... TN,OX,O2,N2,HE,H perturbations from GWs, Nth, East, Ver, Bmer winds
      REAL GWDENTN(9),GWIND(3),GW_BMER   !.. GW added 2023-02-06

     
      REAL UMZ(2,2),AP_TO_KP(401)
      REAL EFAC,EFACN,EFACS,ECDT,ECHROM,ECORON   !.. Solar eclipse factors and time step
      REAL TWISAV                                !.. Save twilight time step
      DOUBLE PRECISION FYLIM,HPLOSS,ZTRAN        !.. For calculating limiting flux
      DIMENSION S(50000),D(19),T(2),FD(9),COSDIP(401),Q(401)
     >   ,V(2),COOL(22)
      COMMON/V91/RVBT(401)
      COMMON/PALT/ZLO,ZHI,JWR,JTI
      COMMON/RFAC/RATFAC(99)
      COMMON/RK/VC(401)
      !----- common to hold the total ion production rates
      REAL SUMION,SUMEXC
      COMMON/BPRTOT/SUMION(3,12,401),SUMEXC(3,12,401)
      !.. N= [O+], [H+], sum minor, [He+] densities, Ti= H+,O+,e temps 
      !.. F= cty eqn functions
      COMMON/FJS/N(4,401),TI(3,401),F(20)
      !.. Ap, solar declination, ephemeris transit, magnetic longitude,
      !.. solar indices, UT in secs, magnetic latitude
      COMMON/MSIS/AP(7),DEC,ETRAN,BLON,F107,F107A,IDAY,SEC,GL(401)
      !.. O+,H+ velocities, field factor (obsolete), magnetic field, 
      !.. GR=gravity, R=geocentric distance,SL= distance along B
      COMMON/VN/U(2,401),BG(401),BM(401),GR(2,401),R(2,401),SL(401)
      !.. [O],[H],[N2],[O2], O+ photoionization rate, Tn
      COMMON/ND/ON(401),HN(401),N2N(401),O2N(401),PHION(401),TN(401)
      !.. used in minor ions [O+], V(O+), [He], ???, He+ prod, ???
      COMMON/A/OPLS(401),VOPLS(401),HE(401),UE(2,401),HEPRD(401),OP(401)
      !.. Altitude(km), time step, ???, ???, ???, ITF= change solution
      !.. mode, JMAX=# grid points, JMAX1=JMAX-1,# Newton iterations, # ions
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
      !.. NSAVE= saved [O+] and [H+], O+, H+ flux, EHT=thermal e heating rate
      COMMON/SAV/NSAVE(2,401),TISAV(3,401),FY(2,401),UN(401),EHT(3,401)
      !.. solar zenith angle, ??? ??? ??? ???, max # time steps ??? ???
      COMMON/LPS/SZA(401),EPSN,DC,ISPC,I1,I2,IPV,JTIMAX,IPP,ISKP
      !.. Stores vibrational N2 and odd N densities for printing. 
      !.. RCON= effect of vibrational N2 on O+ + N2 reaction
      COMMON/STAW/STR(12,401),RCON(401)
      !.. Factors for multiplying EUV, UV, and MSIS
      COMMON/SOL/UVFAC(59),EUV
      !.. RE= Earth radius, (PLAT,PLON) = geographic location of north
      !.. magnetic pole, Pi, angular velocity of Earth 
      COMMON/GPAR/RE,PLAT,PLON,PI,ANGVEL
      !.. Minor ion and neutral densities O+(2D) O+(2P) N(2D) N(4S), NO 
      COMMON/MIB/OP2D(401),OP2P(401),N2D(401),N4S(401),NNO(401)
      !.. densities, velocities, and masses for minor ions solved by Newton
      COMMON/MINOR/XIONN(9,401),XIONV(9,401),XMASS(9)
      !.. Common for O+ - O collision frequencies
      COMMON/CFACT/DHDU_FAC,COLFAC,CFAC(9),ISCH(9)
      !.. Switches for wind type and equations solved
      COMMON/CONFAC/IWIND,IVIBN2,IODDN,INPLUS,IHEP,ITEMP,IHPOP
      !.. minor neutral species on field line grid
      COMMON/VODDN/VN2D(401),VN4S(401),VNNO(401),VO1D(401),IVERT

      !.. JB= grid point of lower boundary, 3 hour Aps, winds flag, SW=MSIS
      !.. switches, ZPROD = upper boundary for EUV ion production
      DATA JB/1/, AP3HR/-1,33000*0/, IWTIM/-1/,SW/15*1,10*1/,ZPROD/700./
      !.. ITSAV ITAU track elapsed time, ISOL= type of solution, HFAC=
      !.. e heating factor, MAXIT= max Newton iters, JTIDT used to set time
      !.. ZLBDY=lower boundary altitude for solving equations, IRIHED IRI/HWM
      !.. wind switch, max solar zenith angle for ion production
      DATA ITSAV,ITAU, ISOL, HFAC, MAXIT, JTIDT, ZLBDY,  IRIHED, SETCHI
     >  / -10000, 0  ,  0  ,  0.0 ,  17  ,   5   , 120. ,  0   ,  1.7/
      !.. north/south lat and long 250 km above station, Z0=absolute lower
      !.. boudary altitude, SCAL, GRD scale grid B and vertical grid points
      DATA GLONN, GLONF, GLATN, GLATF, Z0, SCAL, GRD / 7*0.0/,APFILE/-1/
      DATA JHFLX/0/,WUNIT8/1/,NCOL_EXB/0/, JJJ/0/, ICAUTN/0/
      DATA IVPERP/-1/,SCALK/-1.0/,PCO_MAX/0.0/, ISYM/0/, SECWF8/-999/
      DATA HMF2N,HMF2S/250.0,250.0/, JK,JC/10,10/, NOCON/0/, JBLATPR/0/
      DATA PCO_LIM/1.03/   !.. minimum allowable L-value
      DATA PCO_WIND/1.08/  !.. minimum allowable L-value for vertical grid
      DATA EFAC,EFACN,EFACS,ECDT/4*1.0/   !.. Eclipse factors
      !.. Unit numbers for writing reaction rate warning
      DATA NUMUNIT/9/,UNITNUMS/3,4,6,9,12,13,16,56,66,13*0/
      DATA FEXISTS/-9/    !.. Used in routine TFILE to detect a file
      !.. Conversion factors from Ap to Kp and vice versa.
      !..Kp = 0   0+   1-   1   1+   2-   2   2+   3-   3   3+   4-   4   4+
      !..ap = 0   2    3    4   5    6    7   9    12   15  18   22   27  32
      !..Kp = 5-   5  5+  6-   6  6+   7-   7   7+    8-   8    8+   9-   9      
      !..ap = 39  48  56  67  80  94  111  132  154  179  207  236  300  400
      DATA AP_TO_KP/1*0,1*0,1*0.3,1*0.7,1*1,1*1.3,1*1.7,2*2,3*2.3,3*2.7,
     >  3*3,4*3.3,5*3.7,5*4,7*4.3,9*4.7,8*5,11*5.3,13*5.7,14*6,17*6.3,
     >  21*6.7,22*7,25*7.3,28*7.7,29*8,64*8.3,100*8.7,1*9/
      DATA SATNSAV/24.0/
      PCO=0.0
      JC=1
      RE=6370  ! Apex Re 6378
      ANGVEL=5.2885E-9
      PI=3.14159
      IWW=0
      CALL FLIPIN(S)      !......... initialize arrays

      !.. $$$$$ UPDATE Version # 2023-01-09 here and in PHEAD & PHEADR $$$$$$
      WRITE(6,676)
 676  FORMAT(/6X,'..... FLIP MODEL with APEX Coords: ...2023-02-17 ...'/
     > /,6X,'*** This version takes 2-3 minutes to set up the Apex',
     > /,6X,' coordinates on a basic PC ***'
     > /,6X,'*** Check bottom of *LOG*.txt file for FLIP updates ***'
     > //,6X,'Consult either Phil Richards (pgrichds@gmail.com)'
     > /,6X,'or Pat Dandenault (Patrick.Dandenault@jhuapl.edu) for help'
     > /,6X,'and the appropriate level of recognition for model use.'
     > //,6X,'A FLIP run has 2 phases, The initial FLIP run gives small'
     > /,6X,'summary files of key parameters such as hmF2 and NmF2 and',
     > /,6X,'one large file. By default these output files have names'
     > /,6X,'like FR*.txt. In the 2nd PRINT phase the large FR*.txt is'
     > /,6X,'read by FLIP to produce several comprehensive output',
     > /,6X,'files with altitude profiles. By default, these files',
     > /,6X,'have names like PR*.txt.'//)
      
      CALL OPEN_FILE() !.. added by Xiang Zhao for PC

      !.. KMON= switch to subroutine for printing. ERR flag keeps
      !.. reading the headers until a proper value is obtained. If
      !.. KMON=-1, no need for large binary file (set NOWR8=1)
      NOWR8=0
 137  READ(5,*,ERR=137) KMON
      IF(KMON.EQ.-1) THEN
        NOWR8=1
        KMON=0
      ENDIF
      IF(KMON.NE.0) CALL MONPRN   !.. branch to printing routine

      WRITE(6,173)
 173  FORMAT(/6X,'--------------- BEGIN FLIP RUN ------------------'/)
      
      !.... Check if there is a file with initial densities and temperatures .......
      CALL TFILE(99,FSTART)  !.. See if a file exists
      IF(FSTART.NE.0) THEN
        ISOL=-1   !.. Signals continuation run
        CALL READ_BQ(SCALK,Q)
        CALL RRSTRT(99,PCO,Z0,SCAL,GRD,GLONN,GLONF,GLATN,GLATF,ISOL,UNN
     >  ,UNC,HMF2N,HMF2S,NMF2N,NMF2S,PLAT,PLON)
      ENDIF

      !....... Input parameters to control run. BLATPR 2020-08-30
      CALL DATRD1(ISOL,ITF,ITFS,I1,I2,IPP,MAXIT,JTIMAX,DTSET,DTMAX
     > ,JTIDT,TWIDT,ZLBDY,IREAD,HPEQ,FPAS,HFAC,IVLPS,TSTOP,JWR,ZLO,ZHI,
     > IWIND,IVIBN2,IODDN,INPLUS,IHEP,HEPRAT,ITEMP,IHPOP,IHPKP,BLATPR)
      IF(ITF.EQ.-9) IWIND=0     !.. use Hedin's wind in quick run
      TWISAV=TWIDT

      CALL RATCHK(NUMUNIT,UNITNUMS) !... check reaction rates standard
             
       !.... input new geophysical parameters
       CALL DATRD2(ISOL,IWIND,AP,BLON,F107,F107A,IDAY,SEC,UVFAC)

      !...... input parameters for setting up a new field line
      CALL DATRD3(ISOL,BLON,SEC,DTSET,F107,JMAX,PCO,Z0,SCAL,GRD
     >     ,L_STAG,L_SHORT,L_PHASE,IDAY)
      JK=20       !.. initial hmF2 point in north
      JC=JMAX-20  !.. initial hmF2 point in north

      SECSTART=SEC

       !... make sure run time is reasonable before reading F10.7
       IF(TSTOP+SEC/3600.0.GT.22*8760) THEN
         TSTOP=8785-SEC/3600
         WRITE(6,792) TSTOP
 792     FORMAT(3X,' ??* WARNING: Total run time must be <= 22 years,'
     >     ,1X,'TSTOP reset to',F10.1,' hours')
       ENDIF

      !.. Not possible to get wind from hmF2 at low or high latitudes
      !.. Need to see if files exist and if they have winds, not hmF2
      !.. 2020-07-23
      IF(PCO.LT.1.10.OR.PCO.GT.13) THEN
         IF(IABS(IWIND).GE.1) CALL TEST_WIND_FILES(IWIND,21,22)
      ENDIF

      !..... reading recently added parameters
      CALL DATRD4(ISOL,TIMEON,TIMEOF,AURALT)

      CLOSE(15)!-- close the file of run control parameters - FLIPRUN.DDD

      !...  For user input of Ap and F10.7. Convert Ap to KP for printing.
      KPREAL=0
      IF(AP(1).GT.400) AP(1)=400
      IF(AP(1).GE.0) KPREAL=AP_TO_KP(NINT(AP(1)+1))

      !... If AP(1) is negative calculate other AP from file for MSIS
      AP(2)=AP(1)
      IF(AP(1).GT.0) AP(7)=-1.0  !.. Flag indicating constant AP, F10.7
      IF(AP(1).LT.0) THEN
        CALL GET_F107AP(IDAY,SEC/3600.0,REAL(TSTOP),AP,SW,F107,F107A,
     >     KPREAL)
        APFILE=1 !.. turn on Ap from file
      ENDIF

      !.. Recalculate UVFAC and flux for variable Ap
      IF(NINT(UVFAC(58)).EQ.-1.OR.NINT(UVFAC(58)).EQ.-3)
     >  CALL FACEUV(F107,F107A,UVFAC,ZFLUX)
      IF(NINT(UVFAC(58)).EQ.-2.OR.NINT(UVFAC(58)).EQ.-4)
     >  CALL FACSEE(F107,F107A,UVFAC,ZFLUX)

      IF(ISOL.EQ.-1) ICAUTN=1
      IF(ISOL.EQ.-1) WRITE(6,820)
 820  FORMAT(/9X,' ...... Continuation of a previous FLIP run ......'/)
      IF(ISOL.EQ.0) WRITE(6,821)
 821  FORMAT(/2X,' *** Creating a new field line. Allow a few hours to
     > settle down ***')

      SECSAV=SEC
      TIMSAV=SEC
      DT=DTSET
      JMAX1=JMAX-1
      NEQ=2*JMAX-4
      IF(Z0.GE.180.AND.IVLPS.GT.0) IVLPS=0
      JQ=(JMAX+1)/2

       !..... Get solar declination and ephemeris transit time
       CALL SOLDEC(IDAY,SEC/3600.0,DEC,ETRAN)

      !..Initialize N2* cooling rate factors
      DO J=1,JMAX
        RVBT(J)=1.0
      ENDDO

      !... See if there is any ExB drift. If so, no oddN nor N2*
      CALL AMBS(1,1.3D4,0.0D0,D,T,CHI,SATEQ,GLON,GLATJ)
      IF(NINT(L_STAG).NE.0) THEN
        !.. Get maximum L at dusk for setting up proper altitude grid
        CALL GET_VPERP(IVPERP,DT,SEC/3600.,SATEQ,TSTOP,PCO,BLON,L_STAG
     >     ,L_SHORT,L_PHASE,VPERPN,NCOL_EXB,UTORLT,PCO_MIN,PCO_MAX,
     >     PCO_LIM)
         !.. warning for low L values
         IF(PCO_MIN.LT.PCO_LIM)  THEN
           WRITE(6,371) PCO_MIN,PCO_LIM
 371       FORMAT(3X,'*** WARNING: ExB drift minimum L-value =',F6.3
     >      ,1X,'is too small. Minimum L-value is reset to',F6.3/)
         ENDIF
         !.. warning for low L values and winds from hmF2
        IF(PCO_MIN.LT.PCO_WIND) THEN
          IF(IWIND.NE.0) WRITE(6,968) PCO_WIND
 968      FORMAT(5X,'CANNOT get wind from hmF2. ExB drift minimum' 
     >    ,1X,'L-value is < ',F6.3,'  using HWM-14 instead'/)
          IWIND=0.0
        ENDIF
      ENDIF

      IF(ISOL.EQ.0) CALL SET_HPEQ(PCO,HPEQ)     !.. reset the equatorial [H+]

      !.. Set up the grid points in the orthogonal q coordinates. If ExB 
      !.. drift, the field line is set up for a higher PCO to ensure
      !.. lower boundary does not rise too much
      PCO_start=PCO
      IF(PCO_MAX.GT.1) PCO_start=PCO_MAX
      IF(SCALK.LT.0) CALL QFIELD(JMAX,PCO_start,RE,Z0,SCAL,SCALK,Q)
      CALL WRITE_Q(1,JMAX,SCALK,Q)              !.. put Qs in P*.dat file
      

      !...............  MAIN TIME STEP BEGINS ....................**

      DO JTI=1,JTIMAX+1
      TWIDT=TWISAV
      !.. set the new time step  ** ++ check dL/dt in VPERP
      CALL SETDT(DT,DTSET,DTMAX,JTIDT,ITAU,Z,JMAX,ITF,ITER,TWIDT)
      IF(NINT(DTSET).EQ.59) DT=DTMAX   !..forces constant time steps from start

      !.. Calculate solar eclipse factor in North then South and reset DT
      CALL ECFUN(SEC,IDAY,45.0,ECDT,ECHROM,ECORON)
      EFACN=ECHROM  !.. 2017-09-08: Use chromospheric emission for photoionization
      IF(EFACN.LT.1) WRITE(6,'(7X,2A,F6.3,A,F6.3)') 
     >  'Eclipse in Northern Hemisphere.',
     >  '  Chromosphere irradiance fraction =',ECHROM,
     >  ';  Corona irradiance fraction =',ECORON
      CALL ECFUN(SEC,IDAY,-45.0,ECDT,ECHROM,ECORON)
      EFACS=ECHROM  !.. 2017-09-08: Use chromospheric emission for photoionization
      IF(EFACS.LT.1) WRITE(6,'(7X,2A,F6.3,A,F6.3)') 
     >  'Eclipse in Southern Hemisphere.',
     >  '  Chromosphere irradiance fraction =',ECHROM,
     >  ';  Corona irradiance fraction =',ECORON
      IF(ECDT.LT.DT.AND.(EFACN.LT.1.0.OR.EFACS.LT.1.0)) THEN
        DT=ECDT
        IF(TWIDT.GT.DT) TWIDT=DT  !.. for printinf data each time step for eclipse
      ENDIF

      !--  time to stop? Check # iterations and elapsed time
      IF(JTI.GE.JTIMAX+1) THEN
         WRITE(6,*) ' ** Completed # of time steps'
          CALL FLIP_MODS
          WRITE(6,676)
          STOP 1
       ELSEIF((SEC-SECSTART)/3600.0.GE.TSTOP) THEN
          WRITE(6,*) ' ** Completed the allotted run time'
          CALL FLIP_MODS
          WRITE(6,676)
          STOP !2
      ENDIF

      CALL ACTUAL_DAY(IDAY,SEC,JDAY,SX)    !.. gets actual day of year

      !.... Get the next AP from the file
      IF(APFILE.GE.0)
     >   CALL GET_F107AP(IDAY,SEC/3600.0,REAL(TSTOP),AP,SW,F107,F107A,
     >     KPREAL)
ccc      WRITE(6,'(22F10.1)') F107,F107A,AP

      !... adjust equatorial density for large Kp, only if from file
      IF(APFILE.GT.0.AND.IHPKP.EQ.1)
     >     CALL HP_KP(KPREAL,PCO,SEC/3600.0,JMAX,Z,N)


      !... Get ExB drift velocities. AMBS called for LT at equator
      CALL AMBS(1,1.3D4,0.0D0,D,T,CHI,SATEQ,GLON,GLATJ)
      IF(NINT(L_STAG).NE.0) THEN
         !... determine the phase of teardrop based on Kp/Ap
         !..L_PHASE=TD_PHASE(SEC/3600.0,AP(2),AP(3))
         CALL GET_VPERP(IVPERP,DT,SEC/3600.,SATEQ,TSTOP,PCO,BLON,L_STAG
     >     ,L_SHORT,L_PHASE,VPERPN,NCOL_EXB,UTORLT,PCO_MIN,PCO_MAX,
     >     PCO_LIM)
         !.. Check to adjust DTMAX and TWIDT near the end of run
         CALL ADJDTMAX(ITAU,DT,TSTOP,TWIDT,DTMAX)
      ENDIF

      !... set up the magnetic field line spatial grid .....
      IF(JTI.EQ.1.OR.L_STAG.NE.0) CALL RFIELD(7,PCO,RE,Z0,SCAL,SCALK,
     >    ZLBDY,Q,VPERPN,VTOT,COSDIP,JB,JR,JNQ,JTI,ZHI,ZLO,JVERT,JVERC)


      !... see if auroral file exists and adjust the vertical grid spacing
      CALL TFILE(14,AUR_FILE)
      IF(AUR_FILE.NE.0.AND.GRD.LT.0.5) GRD=1.0
      IF(AUR_FILE.NE.0.AND.Z0.GT.80.0) Z0=80.0

      !--- vertical altitude grid and initial densities for odd N and N2*
      CALL VERGRD(IVLPS,GRD)

      !--- evaluate the geographic location and magnetic declination for headers
      CALL AMBS(JVERT,Z(JVERT),GL(JVERT),D,T,CHIN,SATN,GLONN,GLATN)
      CALL AMBS(JVERC,Z(JVERC),GL(JVERC),D,T,CHIF,SATF,GLONF,GLATF)

      !....... write headers in the files on first time through
      IF(JTI.EQ.1) THEN
         CALL PHEAD(3,PCO,GLONN,GLONF,GLATN,GLATF,Z(JVERT))
         CALL PHEAD(4,PCO,GLONN,GLONF,GLATN,GLATF,Z(JVERT))
         CALL PHEAD(6,PCO,GLONN,GLONF,GLATN,GLATF,Z(JVERT))
         CALL PHEAD(9,PCO,GLONN,GLONF,GLATN,GLATF,Z(JVERT))
         CALL PHEAD(12,PCO,GLONN,GLONF,GLATN,GLATF,Z(JVERT))
         CALL PHEAD(17,PCO,GLONN,GLONF,GLATN,GLATF,Z(JVERT))
         CALL PHEAD(18,PCO,GLONN,GLONF,GLATN,GLATF,Z(JVERT))
         CALL PHEAD(56,PCO,GLONN,GLONF,GLATN,GLATF,Z(JVERT))
         CALL PHEAD(66,PCO,GLONN,GLONF,GLATN,GLATF,Z(JVERT))

         !-- Print the run control parameters - PCOFM obsolete
         CALL RUNPRN(6,ISKP,Z0,HPEQ,FPAS,HFAC,UVFAC,AUR_FILE,PCO)

         !--  Print info about data in wind file and ExB
         CALL PWINDS(IWIND,SATN,SATF,SEC,TSTOP,COSDIP(JK),ZTOPTE,
     >     NCOL_EXB,IODDN,IVIBN2)
         WRITE(6,116) ISOL,ITF,JMAX,IWIND,MAXIT,IVLPS,JTIMAX
     >   ,NINT(HPEQ),FPAS,HFAC,NINT(DTSET),NINT(DTMAX),NINT(TWIDT)
     >   ,TSTOP,HEPRAT,NINT(Z0),COLFAC,(ETRAN-43200.0)/86400.0
      ENDIF

 116  FORMAT(/1X,'ISOL ITF JMAX IWIND MAXIT IVLPS  JTIMAX '
     > ,3X,'HPEQ',2X,'FPAS',1X,'HFAC',2X,'DTSET',2X,'DTMAX',2X
     > ,'TWIDT',2X,'TSTOP',3X,'HEPRAT',5X,'Z0',2X,'COLFAC Eph_tran'
     > ,/2I4,4I6,I9,I8,2F5.1,3I7,2F9.2,I5,F7.1,F7.2)

      CALL AMBS(JVERT,Z(JVERT),GL(JVERT),D,T,CHIN,SATN,GLONN,GLATN)
      CALL AMBS(JVERC,Z(JVERC),GL(JVERC),D,T,CHIF,SATF,GLONF,GLATF)

      !.... Alternative winds are read from file. Otherwise Hedin's winds
      !.... RATIN and RATIS are used to normalize nmF2 to IRI if IWIND<-1
c      WRITE(6,*) ' @#$ call altwin'
      IF(IWIND.NE.0) 
     >  CALL ALTWIN(DT,Z,N,ISOL,IWIND,TIMSAV,TSTOP,IWTIM,IWW,JMAX
     >  ,COSDIP,UNN,UNC,HMF2N,HMF2S,NMF2N,NMF2S,UN,DNMF2N,DNMF2S
     >  ,DHMF2N,DHMF2S,RATIN,RATIS,TOPTEN,TOPTES,ZTOPTE,APFILE)

      ! When fitting topside Te, no plasmaspheric heating by e*
      IF(ZTOPTE.GT.0.0) FPAS=2.0

      IF(ITF.LE.-9) THEN
         WRITE(6,55) SEC/3600.,JDAY,(AP(IAP),IAP=1,7),NINT(F107)
     >      ,NINT(F107A)
         GO TO 221
      ENDIF

      !..  index of print lower bound for electron heating and cooling
      DO J=1,JMAX/2
         IF(Z(J).LT.ZLO) JZLO=J  
      ENDDO
     
55    FORMAT(' UT=',F7.2,', DAY=',I8,'  AP(1..7)= ',7(1X,F5.1)
     >   ,' F10.7=',I4,' F10.7A=',I4)
      WRITE(6,'(A,I5)') ' @#$ DO 17, JTI=',JTI
      !---------- initialization and neutral atmosphere loop ----------------
      DO 17 J=1,JMAX
         !..... o+ + n2* reaction rate factor .....
         VC(J)=1.0
         IF(RCON(J).LT.1.0) RCON(J)=1.0
         IF(RCON(J).GE.1.0) VC(J)=RCON(J)
         IF(JTI.EQ.1.AND.GL(J)*57.2958.GE.BLATPR) JBLATPR=J  !.. 2020-08-30
         IF(J.LE.JQ.AND.Z(J).LT.250.0) J250=J     ! index for F2 region

         IF(J.EQ.9) WRITE(6,'(A,I5)') ' @#$ CALL AMBS',JTI
         !... Obtain neutral densities -- AMBS calls MSIS ......
         CALL AMBS(J,Z(J),GL(J),D,T,CHI,SAT,GLON,GLATJ)
         HE(J)=D(1)
         ON(J)=D(2)
         HN(J)=D(7)
         N2N(J)=D(3)
         O2N(J)=D(4)
         TN(J)=T(2)
         SZA(J)=CHI
         N4S(J)=D(8)
         !.. Evaluate the magnetic field line cmpt. of neutral wind UN(J)
         !.. from Hedin's model (BCOMPU=component along magnetic meridian)
         IF(J.EQ.9) WRITE(6,'(A,I5)') ' @#$ WIND'
          IF(IWIND.EQ.0) THEN
           CALL HEDWIN(J,JDAY,SX,REAL(Z(J)),GLATJ,GLON,SAT,F107A
     >       ,F107,AP,BCOMPU,UHED)
         !... Get GW data to modify winds. 2023-02-06
         !... GWDENTN = TN,OX,O2,N2,HE,H perturbations from GWs   
         IF(J.EQ.9) WRITE(6,'(A,I5)') ' @#$ GW'
         CALL GRAVITY_WAVE_DATA
     >     (0,J,SEC/3600.0,REAL(Z(J)),GLATN,GWDENTN,GWIND,GW_BMER)
           !.. Gravity wave horiz component GW_BMER added 2023-02-06
           UN(J)= -ABS(COSDIP(J))* (BCOMPU+GW_BMER) * 100.0
           IF(J.EQ.JK) UNN= (BCOMPU+GW_BMER) * 100.0
           IF(J.EQ.JC) UNC=-(BCOMPU+GW_BMER) * 100.0

           !.. Output added 2023-02-09
           IF(J.EQ.JZLO.AND.ABS(GWDENTN(2)).GT.0.0) 
     >      WRITE(6,'(A,A,/2X,8F8.2,1P,9E10.2)') 
     >      '      UT    ALT   HWM_Bmer GW_Bmer   GWNth   GWEst   GWVer'
     >      ,'   GW_TN   GW_OX     GW_O2     GW_N2      GW_H     GW_HE'
     >      ,SEC/3600.0,Z(J),BCOMPU,GW_BMER,GWIND(1),GWIND(2),GWIND(3)
     >      ,(GWDENTN(I),I=1,6)
          ENDIF
     
         !..... Calculate centrifugal force and add to gravity
         !..      IF(JTI.EQ.1) CALL CENFOR(IEXP,BLON,GL(J)*57.296,Z(J),CF)
         !..      IF(JTI.EQ.1.AND.IVLPS.GE.0) GR(2,J)=GR(2,J)+CF

         !..... to set up initial temperature and density profiles
         IF(HPEQ.GT.0.0) CALL PROFIN(ISOL,J,JMAX,Z(J),F107,N,TI,TN
     >       ,SL,GR,ON,HN,HPEQ,HEPRAT,XIONN,IVLPS)

         !..... printing of ambient parameters.................
         IF(NINT(SEC/3600.).NE.-43) GO TO 17
         JNEUT=JNEUT+1
         IF(JNEUT.EQ.1) WRITE(3,112)
 112     FORMAT(/'  J      ALT    GL     SZA     GR    TN     N(O)'
     >     ,5X,'N(H)     N(N2)     N(O2)     N(HE)   WIND')
         IF(J.LT.JMAX/2.AND.Z(J).GE.100.AND.Z(J).LE.350)
     >   WRITE(3,'(I4,F8.0,2F7.2,2F7.0,1P,5E9.1,0P,F9.1)')
     >       J,Z(J),GL(J),SZA(J),GR(2,J),TN(J),ON(J),HN(J)
     >      ,N2N(J),O2N(J),HE(J),UN(J)/100.0
 17     CONTINUE
      ENDDO
C ---------------- end neutral atmosphere loop -------------------
      WRITE(6,'(A,I5)') ' @#$ END DO 17'
      
      IF(HPEQ.GT.0) HPEQ=-HPEQ
      !.. Obtain symmetric ambient conditions - diagnostics
      IF(ISYM.NE.0) CALL SYM(JMAX,HE,EHT,UN,Z,SZA)

      !--- transfer vertical grid odd N to field line grid. If no odd N,
      !---- a flag is set in TRSTR 
      CALL TRSTR(IVLPS)

      !.... evaluate primary production rates from primpr
      DO J=1,JMAX
         IF(J.LE.JMAX/2) EFAC=EFACN
         IF(J.GT.JMAX/2) EFAC=EFACS
         CALL PRIMPR(J,Z,ON,N2N,O2N,HE,SZA,TN,JMAX,SETCHI,IVLPS,EFAC)
      ENDDO
      !..  pe2s=2-stream prog to get electron heating rate;
      DO I=1,3
      DO J=1,JMAX
        EHT(I,J)=0.0
      ENDDO
      IF(CHIN.LE.2.0.OR.CHIF.LE.2.0) CALL PE2S(Z,ON,N2N,O2N,HE,BM
     >   ,BG,JMAX,DH,N,TI,EHT,TN,SZA,FPAS,IVLPS,-1.0E22,ISOL,SL)

      !..... If there is an EUV flux file need to reevaluate primary rates
      IF(FEXISTS.EQ.-9) CALL TFILE(42,FEXISTS)
      IF(FEXISTS.NE.0) THEN
        DO J=1,JMAX
          IF(J.LE.JMAX/2) EFAC=EFACN
          IF(J.GT.JMAX/2) EFAC=EFACS
          CALL PRIMPR(J,Z,ON,N2N,O2N,HE,SZA,TN,JMAX,SETCHI,IVLPS,EFAC)
        ENDDO
      ENDIF

      !.... AURSUB calculates auroral electron production rates
      UTHRS=SEC/3600.
      CALL AURSUB(1,JMAX/2,UTHRS,Z,ON,O2N,N2N,N,TI,EHT,DTMAX,AURALT)

      !-- Sum the EUV, photoelectron, and auroral production rates
      CALL SUMPRD(JMAX)

      !.... adjust the plasmaspheric heating rate
      CALL CAL_EHEAT(ISOL,JQ,JHN,JHC,JHX,JMAX,JK,JC,ZTOPTE
     >   ,Z,TN,TI,PCO,HFAC,TIMEON,TIMEOF,F107,TOPTEN,TOPTES
     >   ,SEC,SATN,SATF,SATEQ,TISAV,N,EHT,HEATF)

      WRITE(6,*) ' @#$ after cal_eheat'
      IF(JTI.GT.2) STOP

      !.. loop to calculate ionization rates
      DO J=1,JMAX 
         PHION(J)=1.0E-22
         N(3,J)=1.0E-22
         DO IX=1,9
            FD(IX)=0.0
         ENDDO
         IF(Z(J).GE.Z0.AND.Z(J).LE.ZPROD) THEN
            !.. CALL cminor to get no+, o+(2d), n2+ & o+(2p) densities,
            CALL CMINOR(0,J,0,IPR,FD,7,HE(J))
            !-- PHION=total O+(4S) prod, including EUV, e*, dissoc of O2 and
            !-- minor ions. BCKPRD = small background production to eliminate 
            !-- instability below 200 km
            BCKPRD=2.0E-10*N2N(J)*EXP(-5.0E-14*N2N(J)*TN(J))
            PHION(J)=SUMION(1,7,J)+SUMION(2,4,J)+SUMION(2,5,J)+FD(9)
            SUMION(3,9,J)=FD(9)   !%%%%%%%%%
            PHION(J)=PHION(J)+BCKPRD
         ENDIF
         !.. Sum minor ions N+, NO+, O2+ for electron density at low altitudes
         !.. Leave out N2+, O+(2P), O+(2D) causes [e] gradient problems
         N(3,J)=XIONN(1,J)+FD(1)+FD(2)+(FD(3)+FD(4)+FD(5))*0.0000
      ENDDO

      !... Obtain symmetric ambient conditions - diagnostics
      IF(ISYM.NE.0) CALL SYM(JMAX,HE,EHT,UN,Z,SZA)

       !----- O+, H+, DENSITIES + electron and ion TEMPERATURES ---
       IF(I1.GT.0)  CALL LOOPS(S,NEQ,N,TI,ITAU,JB,SEC,MAXIT,ISOL
     >    ,JK,JC,IWIND,RATIN,RATIS,HE)

       !... testing for nonconvergence
       IF(ITER.EQ.-1) CALL RUN_ERROR    !.. print ERROR warning in output files
       IF(ITER.EQ.-1) STOP 1        !.. could not find a solution

       IF(I1.GT.0.AND.ITER.LE.-2) GO TO 221  !.. jump to end of time loop

      !.. Obtain symmetric ambient conditions - diagnostics
      IF(ISYM.NE.0) CALL SYM(JMAX,HE,EHT,UN,Z,SZA)

      !....... diagnostic print  from subroutines TFIJ and DFIJ ....
      !-- Moved here 21 April 1992 to ensure continuity equation balance
      IF((REAL(ITAU)/3600.0.GE.TSTOP-0.5).AND.JWR.GT.0) THEN
        WRITE(6,*) ' Writing diagnostic info in FOR03*.DAT files'
        IF(JWR.EQ.2) CALL AVDEN(JMAX1,TI,4,ZLO,ZHI)
        DO J=2,JMAX1
           IF(JWR.EQ.1.AND.Z(J).GT.ZLO.AND.Z(J).LT.ZHI)
     >       CALL TFIJ(J,4,IPR,N,TI,F,JB,JBS,V,HE)
           IF(JWR.EQ.2.AND.Z(J).GT.ZLO.AND.Z(J).LT.ZHI)
     >       CALL DFIJ(J,4,IPR,N,TI,F,JB,JBS,V)
         ENDDO
      ENDIF  !.. end print

       !--------------------------- He+ DENSITIES --------------------
c       IF(IABS(IVLPS).GE.9)
        CALL XION(S,NEQ,N,TI,ITAU,JB,SEC,MAXIT,ISOL,9)

       !---------------------------  N+ DENSITIES --------------------
c       IF(IABS(IVLPS).GE.11)
        CALL XION(S,NEQ,N,TI,ITAU,JB,SEC,MAXIT,ISOL,11)

       !... testing for nonconvergence
       IF(ITER.EQ.-1) CALL RUN_ERROR  !.. print ERROR warning in output files
       IF(ITER.EQ.-1) STOP 1        !.. could not find a solution
       IF(I1.GT.0.AND.ITER.LE.-2) GO TO 221

      !---------- N(2D), N(4S), NO, and N2* DENSITIES -----------
      !----- The odd N and N2* time step is never less than TWIDT
      VDT=SEC-SECSAV
      IF(ISOL.EQ.-1.AND.VDT.LT.TWIDT.AND.ITF.EQ.1) GO TO 277
      IF(VDT.LE.0) VDT=DT
      SECSAV=SEC

      IF(IODDN.GE.1.OR.IVIBN2.GE.1)
     >   CALL VERCON(J250,NOCON,S,JMAX,999,VDT,GL(JVERT),GRD,ISOL,IVLPS
     >   ,IVIBN2,HMF2N,HMF2S,SATN,SATF,CHIN,CHIF)

      !... Resetting densities to previous values due to odd N or N2*
      !... non-convergence
      IF(NOCON.EQ.0)  THEN
         CALL SAVDEN(1)       !... save good densities
      ELSE
         CALL RECDEN(ITER)    !... recall previous densities
         GO TO 221
      ENDIF

      IF(ISOL.EQ.0) ITF=ITFS
      ISOL=-1

 277   CONTINUE

      ! ******************** STORE RESULTS *********************
      !....... write a new file at each time step ........
      CALL WRITE_BQ(1,JMAX,SCALK,Q)
      CALL WRSTRT(PCO,Z0,SCAL,GRD,GLONN,GLONF,GLATN,GLATF,UNN,UNC
     >  ,HMF2N,HMF2S,NMF2N,NMF2S,PLAT,PLON)

      !-- Printing summary densities and temperatures in ASCII files
      !-- First, write warning about initial conditions then data in file
      IF(ICAUTN.EQ.0.AND.REAL(ITAU)/3600.0.GT.12) CALL PRCAUT(ICAUTN)
      !.. Get the winds at altitude ZLO in both hemispheres
      CALL HEDWIN(JZLO,JDAY,SX,REAL(ZLO),GLATN,GLONN,SAT,F107A,F107,AP,
     >  BCOMPU,UHED)                 !.. North wind
        UMZ(1,1)=UHED(1)
        UMZ(1,2)=UHED(2)
        JZLOC=JMAX+1-JZLO
      CALL HEDWIN(JZLOC,JDAY,SX,REAL(ZLO),GLATF,GLONF,SAT,F107A,F107,AP,
     >  BCOMPU,UHED)                 !.. South wind
        UMZ(2,1)=UHED(1)
        UMZ(2,2)=UHED(2)
      CALL DIAGPR(ISOL,JDAY,JB,JK,JC,JQ,JNQ,JFQ,HMF2N,HMF2S,UMZ,NMF2N
     >  ,NMF2S,SATN,SATF,CHIN,CHIF,UNN,UNC,VTOT,PCO,VPERPN,COSDIP)

      !-- Output for comparison with IRI model 20 APR 1992
      CALL EDEN170(JMAX,Z,N,J170,N170KMN,N170KMS)  !.. get [e] at 170 km
      CALL IRIPRN(JVERT,JVERC,IRIHED,SEC,SATN,CHIN,NMF2N
     >    ,DNMF2N,HMF2N,DHMF2N,SATF,CHIF,NMF2S,DNMF2S,HMF2S,DHMF2S
     >    ,IDAY,GLATN,GLONN,GLATF,GLONF,F107A,N170KMN,N170KMS)

      !... output transition altitudes 24/2/2000
      CALL TRAN_ALT(N,JMAX,Z,ON,HN,SATN,SATF,UTHRS)

      !-- Output production frequencies
      DO J=1,JMAX/2
        IF(Z(J).LT.600) JFREQN=J
      ENDDO
      JFREQS=JMAX+1-JFREQN
      CALL EUVPRN(JK,JC,JFREQN,JFREQS,SEC,SATN,CHIN,SATF,CHIF,IDAY
     >   ,F107,F107A,TI,Z)
      CALL CAL_STATS(SEC,DT,HMF2N,DHMF2N,HMF2S,DHMF2S)

      !.. Calculate Richards and Torr limiting H+ flux 
      CALL FLXLIMIT(JR,JQ,JK,Z,N,TI,BM,GR,ON,HN,TN,SL,JZ,FYLIM,HPLOSS,
     >   SH_OPLUS,ZTRAN)
      FLIPFLX=N(2,JNQ)*U(2,JNQ)*BM(1)/BM(JNQ)

       !.. Find the peak H+ density to use as transition altitude
      JHP=JQ   !.. begin at hmF2 and go up the field line
      DO J=JQ-1,JK,-1
        !.. Check for peak density altitude
        IF(N(2,J).GT.N(2,J-1).AND.N(2,J).GT.N(2,J+1))
     >    JHP=J
      ENDDO
      
         IF(MOD(JTI-1,18).EQ.0) WRITE(6,317)
 317     FORMAT(/10X,'[------Northern hemisphere-------]',5X,
     >   '[------Southern hemisphere------]',/3X,'  UT    LT  SZA'
     >   ,3X,'NmF2   HmF2  N2* Lat Lon     LT  SZA   NmF2   HmF2'
     >   ,2X,'N2* Lat Lon  Kp   F107   JTI  DelT ITR   L-Sh   YYYYddd')
 117     FORMAT(F8.2,F6.2,I4,1P,E9.2,I5,0P,F5.2,2I4,2X,F6.2,I4,1P
     >   ,E9.2,I5,0P,F5.2,2I4,F5.1,I5,I8,I5,I4,F9.3,I8)
         WRITE(6,117) (SEC)/3600.,SATN,NINT(57.3*CHIN),NMF2N
     >   ,NINT(HMF2N),RCON(JK),NINT(GLATN),NINT(GLONN),SATF
     >   ,NINT(57.3*CHIF),NMF2S,NINT(HMF2S),RCON(JC),NINT(GLATF)
     >   ,NINT(GLONF),KPREAL,NINT(F107),JTI,NINT(DT),ITER,PCO,JDAY

       !.. Special output above an equatorial station 2020-08-30
       !.. Only output if BLATPR is not exactly zero
      IF(JBLATPR.GE.1.AND.JBLATPR.LE.JMAX.AND.ABS(BLATPR).GT.0.001) THEN
        IF(JTI.EQ.1)
     >     CALL PHEAD(43,PCO,GLONN,GLONF,GLATN,GLATF,Z(JVERT))
        CALL AMBS(JBLATPR,Z(JBLATPR),GL(JBLATPR),D,T,CHIX,SATX,
     >       GLONX,GLATX)
        IF(JTI.EQ.1) WRITE(43,'(/5X,A,A,2F5.1,A,F5.1,F6.1,A)') 
     >   ' Ionosphere quantities located on field line above',
     >   ' (MagLat. MagLon.) = (', BLATPR,BLON,'); [GeoLat. Geolon] = ['
     >   ,GLATX,GLONX,']'
        IF(JTI.EQ.1) WRITE(43,'(11X,A,/3X,A,7X,A)') 
     >   'Ne = electron density; min+ = O2+ + NO+ + N2+ + N+'
     >  ,'UT       LT      SZA      ALT     Tn      Ti      Te       Ne'
     >  ,'O+       H+      min+     He+       OX      O2       N2'

        WRITE(43,'(3F8.2,F10.2,3F8.2,1P,9E9.2)') 
     >    SEC/3600.,SATX,57.3*CHIX
     >   ,Z(JBLATPR),TN(JBLATPR),TI(1,JBLATPR),TI(3,JBLATPR)
     >  ,N(1,JBLATPR)+N(2,JBLATPR)+N(3,JBLATPR)+N(4,JBLATPR)
     >  ,N(1,JBLATPR),N(2,JBLATPR),N(3,JBLATPR),N(4,JBLATPR)
     >  ,ON(JBLATPR),O2N(JBLATPR),N2N(JBLATPR)
      ENDIF !.. End special output 2020-08-30
      
      !.... Write topside heat flow and Brace and Theis Te and write
      IF(JHFLX.EQ.0) THEN
         !..JJ=JK    !.. altitude for cooling rates. Can use JC also
         JJ=JZLO  !.. altitude for cooling rates. Can use JK or JC also
         WRITE(56,893) NINT(Z(JHX)),NINT(Z(JJ))
 893     FORMAT(/3X,'Hfac= heating factor, heat_eq= heating rate at'
     >    ,1X,'equator, FLUX= heat flux at',I5,' km (eV cm-2 s-1)'
     >    ,3X,'other heating and cooling rates are at',I7,'km'/
     >    ,3X,'The last 5 colums are for the Richards and Torr 1985'
     >    ,1X,'North limiting H+ flux, HPLOSS = H+ loss above Z0'/)
         IF(NINT(ZTOPTE).GT.200) THEN
           WRITE(56,894) NINT(Z(JHN)),NINT(Z(JHN))
 894       FORMAT(' * Te data constraint imposed at ',I5,' km'
     >           ,1X,' and Tcalc is also at',I5,' km')
      ELSE IF(NINT(ZTOPTE).EQ.1) THEN
           WRITE(56,895)
 895       FORMAT(' * Te data constraint imposed at hmF2'
     >           ,1X,' and Tcalc is also at hmF2')
      ELSE
           WRITE(56,896) NINT(Z(JHX))
 896       FORMAT(' * NO Te data constraint, Tcalc is at',I5,' km')
         ENDIF
      ENDIF    !... end header info


      IF(SATN-SATNSAV.LT.0)          WRITE(56,897)
 897     FORMAT(2X,'  UT      LT     hmF2    SZA Tmeas  Tcalc Ti_cal'
     >    ,2X,'Hfac    FLUX    heat_eq  Tn_hm Ti_hm   e_heat'
     >    ,2X,'ion_cool  N2_cool  O_cool   e_N2D      Ne     H_C'
     >    ,11X,'Z0       ZTRAN   PhiLIM    PhiFLIP   HPLOSS    PhiRAT'
     >    ,5X,'Oplus    Hy_Z0    SH_OPLUS ZHpeak   Tdiff')
      SATNSAV=SATN

          IF(NINT(ZTOPTE).LT.1) JHN=JHX     !... no Te constraint
          JHFLX=1
          HFLUX1=7.7E+5*TI(3,JHX)**2.5*(TI(3,JHX)-TI(3,JHX-1))/
     >      (SL(JHX)-SL(JHX-1))
          !...  DNE=N(1,JHN)+N(2,JHN)+N(3,JHN)+N(4,JHN)
          !... CALL BRACE(Z(JHN),DNE,TI(3,JHN),TB,TEOTB)

          CALL ECRATS(0,JJ,Z(JJ),TI(3,JJ),TI(2,JJ),TN(JJ),N(1,JJ)
     >      ,N(2,JJ),N(3,JJ),ON(JJ),N2N(JJ),O2N(JJ),HN(JJ),HE(JJ)
     >      ,EHT(3,JJ),TLSS,FOL,COOL)
            CALL CMINOR(0,JJ,0,IPR,FD,7,HE(JJ))

          WRITE(56,218) SEC/3600,SATN,Z(JK),NINT(57.3*CHIN)
     >      ,NINT(TOPTEN),NINT(TI(3,JHN)),NINT(TI(1,JHN)),HEATF
     >      ,HFLUX1,EHT(3,JQ),NINT(TN(JJ)),NINT(TI(1,JJ))
     >      ,(COOL(K),K=1,4),FD(8)*2.4,N(1,JJ)+N(2,JJ)+N(3,JJ)
     >      ,COOL(1)+FD(8)*2.4-COOL(5),Z(JZ),ZTRAN,FYLIM
     >      ,FLIPFLX,HPLOSS,FLIPFLX/FYLIM,N(1,JZ),HN(JZ),SH_OPLUS/1E5
     >      ,Z(JHP),(TOPTEN-TI(3,JHN))

 218  FORMAT(3F8.2,4I6,F8.2,1P,E10.2,1P,E10.2,2I6,1P,E10.2,5E9.2,E10.2,
     > 0P,2F10.1,1P,7E10.2,0P,2F8.0)

      !... for printing ion-neutral cooling. not ready yet
      !..        CALL ION_NEUTRAL(0,N(1,JK),N(2,JK),ON(JK),N2N(JK),O2N(JK)
      !..     >   ,HN(JK),HE(JK),TN(JK),TI(2,JK),Z(JK),COOL)

      !.... set a switch to prevent writing into unit 8 until midnight if
      !.... setting up new field line
      IF(ISOL.EQ.0) WUNIT8=0
      !...      IF(ISOL.NE.0.AND.SEC.GE.0.0) WUNIT8=1

      !--- CALL subroutine to write main data on unit 8 at intervals TWIDT
      !--- Data is not written until run has settled down
      IF(ITF.NE.1) GO TO 221
      IF(WUNIT8.NE.0.OR.SATN.LT.2.OR.SATF.LE.2)  THEN
        IF(SEC-SECWF8.GE.TWIDT.AND.NOWR8.EQ.0) THEN
          CALL WF8DAT(IVLPS,JQ,PCO,Z0,SCAL,GRD,GLONN,GLONF,GLATN,GLATF,
     >           ZLBDY,FPAS,UV,HFAC,TIMEON,TIMEOF,AURALT,PLAT,PLON)
          CALL WF8D2(JK,JC,UNN,UNC,IVLPS,JQ,UE)
          SECWF8=SEC
        ENDIF
        WUNIT8=1  !... switches on print from here on
      ENDIF
      ! ******************** END STORE RESULTS *************************

      !... print densities and temperatures from subroutine loops....
      IF(JTIMAX.NE.JTI) GO TO 221
      IF(ITF.EQ.3.OR.ITER.LE.0.OR.JTIMAX.EQ.1) ITSAV=-1000000
      IF(ITAU-ITSAV.LT.3600) GO TO 221
      ITSAV=ITAU
      IF(IPP.EQ.0.OR.JWR.EQ.0) GO TO 221
      WRITE(3,859)
      WRITE(9,859)
 859  FORMAT(/'  J ',3X,'ALT',3X,'TE',4X,'TI',3X,'(O+)',5X,'(H+)',6X
     > ,'VO+',6X,'VH+',5X,'(He+)',5X,'VHe+',5X,'EHT    (N+)     VN+')
      DO J=1,JQ
      JU=JMAX+1-J
      IF(Z(J).LT.ZLO.OR.Z(J).GT.ZHI) GO TO 854
      IF(J.EQ.1.OR.J.EQ.JQ) GO TO 852
      IF(MOD(J,IPP).EQ.0) GO TO 852
      GO TO 854
 852  CONTINUE

      NE=N(1,J)+N(2,J)+N(3,J)
      CALL BRACE(Z(J),NE,TI(3,J),TB,TEOTB)
      WRITE(3,888) J,NINT(Z(J)),NINT(TI(3,J)),NINT(TI(1,J))
     > ,N(1,J),N(2,J),U(1,J),U(2,J),N(4,J),UE(1,J),EHT(3,J)
     > ,XIONN(1,J),XIONV(1,J)
      NE=N(1,JU)+N(2,JU)+N(3,JU)

      CALL BRACE(Z(JU),NE,TI(3,JU),TB,TEOTB)
      WRITE(9,888) JU,NINT(Z(JU)),NINT(TI(3,JU)),NINT(TI(1,JU))
     > ,N(1,JU),N(2,JU),U(1,JU),U(2,JU),N(4,JU),UE(1,JU),EHT(3,JU)
     > ,XIONN(1,JU),XIONV(1,JU)
 854     CONTINUE
      ENDDO
      ENDDO
 221  CONTINUE

 64   FORMAT(2I6,3D12.4,/11F10.2)
 65   FORMAT(16A8)
 66   FORMAT(32A4)
 68   FORMAT(22I6)
 111  FORMAT(I4,F8.0,2F7.2,2F7.0,1P,22E9.1)
 167  FORMAT(10E10.2)
 888  FORMAT(I4,I6,I6,I5,1P,2E9.2,1P,5E9.1,1P,12E9.1)
      STOP
      END

C::::::::::::::::::::::: CAL_EHEAT ::::::::::::::::::::::::::::
      !.. Adjusting plasmaspheric heating rate.
      SUBROUTINE CAL_EHEAT(ISOL,JQ,JHN,JHC,JHX,JMAX,JK,JC,ZTOPTE
     >   ,Z,TN,TI,PCO,HFAC,TIMEON,TIMEOF,F107,TOPTEN,TOPTES,SEC
     >   ,SATN,SATS,SATEQ,TISAV,N,EHT,HEATF)
      IMPLICIT NONE
      INTEGER ISOL,JQ,JHN,JHC,JHX,JMAX,IL,JC,JK,J,IH
      REAL ZTOPTE,ZHFLX,TIMEON,TIMEOF,F107,SATN,SATS,PROP,SATEQ
     >   ,TMAX,SIGT,LMAX,SIGL,FAC_RING,RING_CURR,AMP_RING
     >   ,RC_TIME_CON,SEC
      DOUBLE PRECISION Z(401),TN(401),TI(3,401),HFAC,EHT(3,401),PCO,
     >  TOPTEN,TOPTES,HEATF,TISAV(3,401),N(4,401),HEATFS,TOPTENSV,TETEMP
      !... parameters for the ring current model
      COMMON/RC/AMP_RING,TMAX,SIGT,LMAX,SIGL,RC_TIME_CON,PROP
      DATA ZHFLX/486.0/     !.. altitude of topside heat flux
      DATA HEATFS/0.0/
      
      !.. First determine grid point where Te normalization will be done
      JHN=JMAX/2
      JHC=JMAX/2
      JHX=JMAX/2
      DO 545 IL=2,JMAX/2
        IF(ZTOPTE.LT.(Z(IL+1)+Z(IL))/2.AND.ZTOPTE.GE.(Z(IL)+Z(IL-1))/2)
     >     JHN=IL
        IF(ZHFLX.LT.(Z(IL+1)+Z(IL))/2.AND.ZHFLX.GE.(Z(IL)+Z(IL-1))/2)
     >     JHX=IL
 545  CONTINUE

      JHC=JMAX+1-JHN            !..... Conjugate index
      !... reset if topside Te is at hmF2
      IF(JK.GT.1.AND.NINT(ZTOPTE).EQ.1) JHN=JK
      IF(JK.GT.1.AND.NINT(ZTOPTE).EQ.1) JHC=JC

      !... determine factor for ring current heating of ions
      IF(AMP_RING.GT.0.0.AND.RC_TIME_CON.GT.0.0) THEN
        ZTOPTE=-1
        FAC_RING=AMP_RING*RING_CURR(SATEQ,REAL(PCO),TMAX,SIGT,LMAX,SIGL)
        FAC_RING=FAC_RING*EXP(-SEC/3600.0/RC_TIME_CON)   !.. time decay
        !..WRITE(6,95) NINT(TI(1,JQ)),NINT(TI(3,JQ))
        !.. >     ,NINT(TN(JQ)),FAC_RING,FAC_RING*(N(2,JQ)+N(1,JQ))
c 95     FORMAT(' RC ',3I9,1P,2E11.3) 
        DO J=1,JMAX
          IF(Z(J).GT.2000) THEN
             EHT(1,J)=FAC_RING*(1.-PROP)*(N(2,J)+N(1,J))+EHT(1,J)
             EHT(3,J)=FAC_RING*PROP*(N(2,J)+N(1,J))+EHT(3,J)
          ENDIF
        ENDDO
      ENDIF
      
      IH=0  !.. switch to turn on heating
      !.. before midnight
      IF(HFAC.GT.0.01.AND.SATEQ.GE.TIMEON.AND.SATEQ.LE.TIMEOF) IH=1
      !.. after midnight assuming TIMEOF >= 24
      IF(HFAC.GT.0.01.AND.TIMEOF.GE.24.AND.SATEQ+24.LE.TIMEOF) IH=1
      IF(IH.EQ.1) THEN
        !... Unknown heat source. HFAC=1 -> HF=1.0E9 eV cm-3s-1 OR 
        !... if <0 a straight multiplicative factor at all altitudes
         DO J=1,JMAX
           IF(Z(J).GT.1000)  EHT(3,J)=2.5*HFAC/PCO**4 + EHT(3,J)
         ENDDO
         ZTOPTE=-1  !... signal to return
      ELSE
        !.... If HFAC < 0.0 use constant scaling
        IF(HFAC.LT.-0.01) THEN
          DO J=1,JMAX
           IF(ZTOPTE.GT.200.AND.Z(J).LT.ZTOPTE) EHT(3,J)=-HFAC*EHT(3,J)
           IF(ZTOPTE.LE.200) EHT(3,J)=-HFAC*EHT(3,J)
          ENDDO
           ZTOPTE=-1  !... signal to return
        ENDIF
      ENDIF

      IF(NINT(ZTOPTE).LE.0) RETURN
      IF(ISOL.NE.-1) RETURN  !.. only for regular run

      !.. Increased plasmaspheric heating. HEATF is used when reproducing
      !.. Te input from a file using heat flux algorithm
      !.. 2020-09-27: Added HEATFS and TOPTENSV. Also since the algorithm 
      !.. needs a Te difference, a small percentage is added to TOPTEN
      IF(ZTOPTE.GT.0) THEN
         !... HEATFS, TETEMP, and TOPTENSV smooth out oscillations in HEATF
         HEATFS=HEATF
         IF(TOPTEN.GT.TN(JHN)) 
     >     TETEMP=0.333*TOPTENSV+0.666*TOPTEN+0.04*ABS(TOPTEN-TN(JHN))
         IF(TOPTEN.GT.TN(JK)) HEATF=5E-8*(TETEMP**2.5-TI(3,JHN)**2.5)
         IF(TOPTES.GT.TN(JC)) HEATF=5E-8*(TOPTES**2.5-TI(3,JHC)**2.5)

         !.. make sure heating factor is reasonable. It can be negative in 
         !.. case a previous HEATF was too big.
         IF(HEATF.GT.15.0) HEATF=15.0
         IF(HEATF.LT.-2.0) HEATF=-2.0  

         HEATF=(HEATF+HEATFS)/2.0      !.. Smooth heat flux - not needed
         !.. If MSIS Tn and TOPTEN don't match could get weird profiles 
         !.. and non-convergence
         IF(HEATF.LT.0.0.AND.TI(3,JQ).LT.TETEMP+100) HEATF=0

         TOPTENSV=TETEMP   !.. save for smoothing
         
         !... Now adjust electron EHT(3,J) and ion heating rates
         DO J=1,JMAX
           IF(Z(J).GT.ZTOPTE+500) THEN
              EHT(1,J)=9.0*(1.-PROP)*HEATF/PCO**4 + EHT(1,J)
              EHT(3,J)=9.0*PROP*HEATF/PCO**4 + EHT(3,J)
           ENDIF
         ENDDO
      ENDIF
      
      RETURN
      END
C::::::::::::::::::::::::::: RING_CURR ::::::::::
C... This function fits the ring current using two Guassian distributions
      REAL FUNCTION RING_CURR(LT,       !.. local time
     >                         L,       !.. L-value
     >                      TMAX,       !.. time of max ring current
     >                      SIGT,       !.. width of time Guassian
     >                      LMAX,       !.. L-value of max ring current
     >                      SIGL)       !.. width of L Guassian
      IMPLICIT NONE
      REAL LT,L,TMAX,SIGT,LMAX,SIGL,DELTIME
      !... First term is for time. Second term is for L
      DELTIME=ABS(LT-TMAX)
      IF(ABS(LT-TMAX).GT.ABS(LT+24-TMAX)) DELTIME=ABS(LT+24-TMAX)
      RING_CURR=EXP(-(DELTIME/SIGT)**2) * EXP(-((L-LMAX)/SIGL)**2)
      RETURN
      END
C:::::::::::::::::::::: BLOCK DATA :::::::::::::::::::::::::::::
C..... This BLOCK DATA is used to initialize the arrays in the FLIP
C..... model. This is necessary for the IBM machine
      BLOCK DATA BLO3
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      INTEGER ION
      REAL SEC,DEC,ETRAN,F107,F107A,AP,BLON,EUV,UVFAC,ZPAS,PASK,FPAS
      COMMON/PALT/ZLO,ZHI,JWR,JTI
      COMMON/MSIS/AP(7),DEC,ETRAN,BLON,F107,F107A,IDAY,SEC,GL(401)
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
      COMMON/LPS/SZA(401),EPSN,DC,ISPC,I1,I2,IPV,JTIMAX,IPP,ISKP
      COMMON/SOL/UVFAC(59),EUV
      COMMON/PANGS/ZPAS,PASK,IPAS,IPASC,FPAS
      DATA ZPAS /2000.0/, ZLO/0.0D0/,ZHI/0.0D0/,JWR/0/,JTI/0/
      DATA AP/7*0.0/,DEC/0.0/,ETRAN/0.0/,BLON/0.0/,F107/0.0/,F107A/0.0/
      DATA IDAY/0/,SEC/0.0/, DT/0.0D0/,DH/0.0D0/, ITF/0/,JMAX1/0/
      DATA IPV/0/,ITER/0/, JTIMAX/0/, EUV/0.0/, ISKP/1/, ISPC/5/
      DATA EPSN,  DC,I1,I2,IPP,    THF,EPS  ,ION, TF
     >   /  .999,.5D0, 1, 2, 4 , 1.0D0 ,1E-3, 2 ,1.0 /
      END

C::::::::::::::::::::::::::::::::: ACTUAL_DAY:::::::::::::::::::::::::
C...... Since IDAY does not change and UT continually increases, need to
C...... evaluate actual day (ID) and UT (SX) for MSIS, HWM, and IRI.
C...... Modified by P. Richards in November 2008 to allow for long runs
      SUBROUTINE ACTUAL_DAY(IDAY,    !.... FLIP start day YYYYDDD
     >                       SEC,    !.... Current UT in secs (continuous)
     >                        ID,    !.... actual day YYYYDDD (output) 
     >                        SX)    !.... UT (0-24) in secs (output)
      
      IMPLICIT NONE
      INTEGER SDIM                  !.. Dimension for array
      PARAMETER(SDIM=11)
      INTEGER I,IWR                 !.. Loop control, IWR= write control 
      INTEGER IDAY,ID               !.. FLIP start YYYYDDD, ID= current YYYYDDD
      INTEGER YSTART,DSTART         !.. Start year YYYY and start day DDD
      INTEGER NYRS,NDAYS,NSECS(SDIM)  !.. # of years, days, and seconds in these years
      INTEGER IDAYS                 !.. # Number of days in run so far
      REAL SEC,SX                   !.. UT in secs, actual UT in secs (output)
      DATA IWR/0/

      !.. Check the run is not too long
      IF(SEC.GT.SDIM*3.16E+7) THEN
         WRITE(6,55) SDIM
         CALL RUN_ERROR    !.. print ERROR warning in output files
         STOP
      ENDIF
 55   FORMAT(2X,'** ERROR: FLIP cannot do more than',I3,' year run **')

      YSTART=IDAY/1000               !.. Determine the start year YYYY
      DSTART=MOD(IDAY,1000)          !.. Determine the start day DDD

      !.. Determine the number of days and seconds simulated so far 
      !.. in order to determine the current year and the UT in that year
      NDAYS=0
      DO I=1,SDIM
        IF(MOD(YSTART+I-1,4).EQ.0) THEN     !.. leap year
           NDAYS=NDAYS+366
        ELSE
           NDAYS=NDAYS+365                  !.. normal year
        ENDIF

        !.. Time elapsed since the start day in seconds
        NSECS(I)=(NDAYS-DSTART+1)*24*3600

        !.. Jump out of loop when the time exceeds the current FLIP run time
        IF(NSECS(I).GT.INT(SEC)) THEN
           NYRS=I-1
           GOTO 10
        ENDIF
      ENDDO

 10   CONTINUE

      !.. Determine year and UT in the first year and then later years
      IF(NYRS.LE.0) THEN
        IDAYS=INT(SEC)/3600/24
        ID=(YSTART)*1000+DSTART+IDAYS
      ELSE                                   !.. after first year
        IDAYS=(INT(SEC)-NSECS(NYRS))/3600/24
        ID=(YSTART+NYRS)*1000+IDAYS+1
      ENDIF

      SX=AMOD(SEC,8.64E+4)  !.. UT 

      IF(MOD(ID/1000,4).EQ.0) THEN     !.. leap year
        IF(IDAYS.GT.366) THEN
          WRITE(6,*) '  ** ERROR > 366 days in leap year'
          CALL RUN_ERROR    !.. print ERROR warning in output files
        ENDIF
      ELSE
        IF(IDAYS.GT.365) THEN
          WRITE(6,*) '  ** ERROR > 365 days in normal year'
          CALL RUN_ERROR    !.. print ERROR warning in output files
        ENDIF
      ENDIF

      IWR=0  !..IWR+1  !.. IWR=0 switches off DEBUG write statements
      IF(IWR.EQ.1) WRITE(6,191) 
 191  FORMAT(5X,'IDAY',8X,'ID    NYRS IDAYS     SX',9X,'SEC',8X,'Days')
      IF(IWR.GT.0) WRITE(6,'(2I11,2I5,9I11)') IDAY,ID,NYRS,IDAYS,
     >    INT(SX),INT(SEC),NDAYS,NSECS(NYRS)

      RETURN
      END
C::::::::::::::::::::::::::::: TRAN_ALT :::::::::::::::::::::::::::::
C... This function calculates the transition altitudes from O+ to light
C... ions. Modified 2007-10-18 by P. Richards to do both hemispheres
      SUBROUTINE TRAN_ALT(N,    !.. ion densities
     >                 JMAX,    !.. number of grid points
     >                  ALT,    !.. altitude array
     >                   ON,    !.. oxygen array
     >                   HN,    !.. hydrogen array
     >                 SATN,    !.. Local Time in the Northern hemisphere
     >                 SATF,    !.. Local Time in the Southern hemisphere
     >                 UTHRS)   !.. Universal time
      IMPLICIT NONE
      INTEGER JTOPHPN,JTOPHPC,JTOPLIN,JTOPLIC  !.. transition indices
      INTEGER JMOLN,JMOLC                      !.. molecular transition indices
      INTEGER JMAX                        !.. number of grid points
      INTEGER J,JU,IWR                    !.. loop control variables
      REAL SATN,SATF,UTHRS                !.. Defined above
      REAL RAT1,RAT2                      !.. For calculating ratios
      DOUBLE PRECISION N(4,401),ALT(401)  !.. defined above
      DOUBLE PRECISION ON(401),HN(401)    !.. defined above
      
      !.. Determine the indices where the ion transitions take place
      DO 20 J=1,JMAX/2
        JU=JMAX+1-J                         !.. conjugate index
        IF(N(1,J).GT.N(2,J)) JTOPHPN=J
        IF(N(1,JU).GT.N(2,JU)) JTOPHPC=JU
        IF(N(1,J).GT.N(2,J)+N(4,J)) JTOPLIN=J
        IF(N(1,JU).GT.N(2,JU)+N(4,JU)) JTOPLIC=JU
        IF(ALT(J).LT.500.AND.N(1,J).LT.N(3,J)) JMOLN=J
        IF(ALT(JU).LT.500.AND.N(1,JU).LT.N(3,JU)) JMOLC=JU
 20   CONTINUE

      !------------ Northern Hemisphere --------------------
      !.. See if the next point will give a O+/H+ ratio closer to 1
      RAT1=N(1,JTOPHPN)/N(2,JTOPHPN) 
      RAT2=N(1,JTOPHPN+1)/N(2,JTOPHPN+1)
      IF(ABS(1-RAT1) > ABS(1-RAT2)) JTOPHPN=JTOPHPN+1

      !.. See if the next point will give a O+/(H+ + He+) ratio closer to 1
      RAT1=N(1,JTOPLIN)/(N(2,JTOPLIN)+N(4,JTOPLIN)) 
      RAT2=N(1,JTOPLIN+1)/(N(2,JTOPLIN+1)+N(4,JTOPLIN+1))
      IF(ABS(1-RAT1) > ABS(1-RAT2)) JTOPLIN=JTOPLIN+1

      !------------ Southern Hemisphere --------------------
      !.. See if the next point will give a O+/H+ ratio closer to 1
      RAT1=N(1,JTOPHPC)/N(2,JTOPHPC) 
      RAT2=N(1,JTOPHPC+1)/N(2,JTOPHPC+1)
      IF(ABS(1-RAT1) > ABS(1-RAT2)) JTOPHPC=JTOPHPC+1

      !.. See if the next point will give a O+/(H+ + He+) ratio closer to 1
      RAT1=N(1,JTOPLIC)/(N(2,JTOPLIC)+N(4,JTOPLIC)) 
      RAT2=N(1,JTOPLIC+1)/(N(2,JTOPLIC+1)+N(4,JTOPLIC+1))
      IF(ABS(1-RAT1) > ABS(1-RAT2)) JTOPLIC=JTOPLIC+1
     
      DATA IWR/0/  !.. IWR is used to write the header only once
      IWR=IWR+1
      IF(IWR.EQ.1) WRITE(18,95)
 95   FORMAT(/'   ---- approximate transition altitudes for O+ to H+'
     > ,1X,'and light ions (H+ + He+) and ion densities at those'
     > ,1X,'altitudes in both north and south hemispheres'
     > /,11X,'North Hemisphere O+ to H+ transition   North O+ to' 
     > ,1X,'(H+ + He+) transition',3X,'South Hemisphere O+ to H+'
     > ,1X,'transition    South O+ to (H+ + He+) transition'   
     > /,5X,'UT    LT     ALT      O+      H+      He+      ALT'
     > ,5X,'O+      H+       He+      LT     ALT     O+      H+'
     > ,7X,'He+      ALT     O+      H+       He+      OpMN   OpMS'
     > ,2X,'H/O_Nth  H/O_Sth') 

      WRITE(18,90) UTHRS,SATN,NINT(ALT(JTOPHPN)),N(1,JTOPHPN)
     > ,N(2,JTOPHPN),N(4,JTOPHPN),NINT(ALT(JTOPLIN)),N(1,JTOPLIN)
     > ,N(2,JTOPLIN),N(4,JTOPLIN)
     > ,SATF,NINT(ALT(JTOPHPC)),N(1,JTOPHPC),N(2,JTOPHPC),N(4,JTOPHPC)
     > ,NINT(ALT(JTOPLIC)),N(1,JTOPLIC),N(2,JTOPLIC),N(4,JTOPLIC),
     >  NINT(ALT(JMOLN)),NINT(ALT(JMOLC)),
     >  HN(JTOPHPN)/ON(JTOPHPN),HN(JTOPHPC)/ON(JTOPHPC)
 90   FORMAT(F8.2,F7.2,I7,1P,3E9.2,I7,1P,3E9.2,0P,F7.2,I7,1P,
     >    3E9.2,I7,1P,3E9.2,2I7,1P,2E9.1)
      RETURN
      END
C:::::::::::::::::: HP_KP ::::::::::::::::::::::::
C... Adjust the H+ density for high Kp.
C... Note that Binsack JGR 1967, page 5321 gives Lpp=6.0-0.6*Kp
C... Alternative from Moldwin et al. AGU Spring 2001 SM21A-11
C... Lpp(day)=4.7-0.31*Kp; Lpp(night)=5.5-0.37*Kp
C... Modified July 2008 to allow variable equatorial density instead of 300
      SUBROUTINE HP_KP(KP,      !... Kp index    
     >                PCO,      !... L-shell
     >              UTHRS,      !... Universal time
     >               JMAX,      !... max index on field line
     >                  Z,      !... altitude array
     >                  N)      !... ion densities
      INTEGER J,JMAX
      REAL UTHRS,ALPHA,KP,HP_DEP
      DOUBLE PRECISION PCO,N(4,401),Z(401)

      !.. equatorial density for a depleted flux tube. Note that plasma may be convected to 
      !.. larger L-shells but then some of it convected back from larger L-shells as activity abates
      !.. HP_DEP=30+1.22E4/PCO**4  !.. Old value replaced in August 2012 
      !.. HP_DEP=500.0-100.0*PCO   !.. gives a more gradual decreas with L
      !.. HP_DEP=240.0/PCO/PCO     !.. before 2022-12-20
      HP_DEP=3000.0/PCO/PCO        !.. Updated 2022-12-20
      IF(HP_DEP.LT.5.0) HP_DEP=5.0   !.. set floor value Chappell 1972, fig.4


      !.. No need to deplete already depleted flux tube
      IF(N(2,(JMAX+1)/2).LT.HP_DEP) RETURN

      ALPHA=HP_DEP/N(2,(JMAX+1)/2)   !.. reduction fac. at equator

      !.. If Kp is high enough, reduce H+ and He+ densities. Note that 24 hrs
      !.. must have elapsed since previous reduction, otherwise too unstable. 
      IF(KP.GE.((6.0-PCO)/0.6).AND.UTHRS-THPEQ.GT.24.) THEN !..Binsack
        DO J=1,JMAX
          N(2,J)=N(2,J)*ALPHA
          N(4,J)=N(4,J)*ALPHA
        ENDDO
        THPEQ=UTHRS
        WRITE(6,661) N(2,(JMAX+1)/2),KP,UTHRS
      ENDIF
 661  FORMAT(' H+ and He+ reduced; equatorial [H+]=',F10.2
     >   ,1X,' for large Kp=',F4.1,' at UT=',F8.2)
      RETURN
      END 
C:::::::::::::::::::::::::::::: SAVDEN ::::::::::::::::::::::::::::::::
C..... This routine saves all the values to be used after non-convergence in
C..... Odd N solution
      SUBROUTINE SAVDEN(IFLAG)
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      COMMON/ODNSAV/ NSAVE(4,401),TISAV(3,401),XISAV(9,401),XVSAV(9,401)
     >  ,USAV(2,401),UESAV(2,401),STRSAV(12,401),RSAV(401)
      COMMON/FJS/N(4,401),TI(3,401),F(20)
      COMMON/VN/U(2,401),BG(401),BM(401),GR(2,401),R(2,401),SL(401)
      COMMON/A/OPLS(401),VOPLS(401),HE(401),UE(2,401),HEPRD(401),OP(401)
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
      COMMON/STAW/STR(12,401),RCON(401)
      COMMON/MINOR/XIONN(9,401),XIONV(9,401),XMASS(9)

      DO J=1,JMAX
        DO I=1,4
           NSAVE(I,J)=N(I,J)
        ENDDO
        DO I=1,3
          TISAV(I,J)=TI(I,J)
        ENDDO
        DO I=1,2
           UESAV(I,J)=UE(I,J)
           USAV(I,J)=U(I,J)
        ENDDO
        DO I=1,12
           STRSAV(I,J)=STR(I,J)
        ENDDO
        RSAV(J)=RCON(J)
        XISAV(1,J)=XIONN(1,J)
        XVSAV(1,J)=XIONV(1,J)
      ENDDO
      RETURN
      END
C:::::::::::::::::::::::::::::: SAVDEN ::::::::::::::::::::::::::::::::
C..... This routine recovers all the values to be used after non-convergence in
C..... Odd N solution
      SUBROUTINE RECDEN(IFLAG)
      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      COMMON/ODNSAV/ NSAVE(4,401),TISAV(3,401),XISAV(9,401),XVSAV(9,401)
     >  ,USAV(2,401),UESAV(2,401),STRSAV(12,401),RSAV(401)
      COMMON/FJS/N(4,401),TI(3,401),F(20)
      COMMON/VN/U(2,401),BG(401),BM(401),GR(2,401),R(2,401),SL(401)
      COMMON/A/OPLS(401),VOPLS(401),HE(401),UE(2,401),HEPRD(401),OP(401)
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
      COMMON/STAW/STR(12,401),RCON(401)
      COMMON/MINOR/XIONN(9,401),XIONV(9,401),XMASS(9)

      IFLAG=-9       !... flag indicates retrieving stored values
      DO J=1,JMAX
        DO I=1,4
           N(I,J)=NSAVE(I,J)
        ENDDO
        DO I=1,3
          TI(I,J)=TISAV(I,J)
        ENDDO
        DO I=1,2
           UE(I,J)=UESAV(I,J)
           U(I,J)=USAV(I,J)
        ENDDO
        DO I=1,12
           STR(I,J)=STRSAV(I,J)
        ENDDO
        RCON(J)=RSAV(J)
        XIONN(1,J)=XISAV(1,J)
        XIONV(1,J)=XVSAV(1,J)
      ENDDO
      RETURN
      END
C::::::::::::::::::::::::: EDEN170 ::::::::::::::::::::
      SUBROUTINE EDEN170(JMAX,     !.. # points on grid
     >                   Z,        !.. Altitude array
     >                   N,        !.. Ion density array
     >                   J170,     !.. index of 170 km alt in Z
     >                   N170KMN,  !.. [e] at 170 km in north
     >                   N170KMS)  !.. [e] at 170 km in south
      IMPLICIT NONE
      INTEGER J,JMAX,J170,J170C
      DOUBLE PRECISION Z(401),N(4,401),N170KMN,N170KMS
      J170=2 
      DO J=2,JMAX/2
        IF(170.LT.(Z(J+1)+Z(J))/2.AND.170.GE.(Z(J)+Z(J-1))/2) J170=J
      ENDDO
      N170KMN=N(1,J170)+N(2,J170)+N(3,J170)+N(4,J170)
      J170C=JMAX-J170+1
      N170KMS=N(1,J170C)+N(2,J170C)+N(3,J170C)+N(4,J170C)
      RETURN
      END
C::::::::::::::::::::::::: EDEN170 ::::::::::::::::::::
C.. This routine calculates the sum of the squares of the differences between the
C.. measured and modeled hmF2
      SUBROUTINE CAL_STATS(SEC,DT,HMF2N,DHMF2N,HMF2S,DHMF2S)
      IMPLICIT NONE
      INTEGER JTI
      REAL SEC
      DOUBLE PRECISION HMF2N,DHMF2N,HMF2S,DHMF2S,SUMSQN,SUMSQS,NUMSQ,DT,
     >  SUMN,SUMS
      DATA JTI/0/

      IF(JTI.GT.-1.1111) RETURN  !...

      JTI=JTI+1
      !.. zero out every 24 hours to see problem data
      IF(JTI.EQ.1.OR.AMOD(SEC,86400.0).LT.DT+1) THEN
       IF(NUMSQ.GT.0) WRITE(6,441) 
     >   NINT(NUMSQ),DSQRT(SUMSQN/NUMSQ),DSQRT(SUMSQS/NUMSQ),
     >   SUMN/NUMSQ,SUMS/NUMSQ
 441    FORMAT(' #=',I7,' SIGMAN=',1PE10.2,' SIGMAS=',1PE10.2,
     >   ' SUMN=',1PE10.2,' SUMS=',1PE10.2)
        SUMSQN=0.0
        SUMSQS=0.0
        NUMSQ=0.0
        SUMN=0.0
        SUMS=0.0
      ENDIF

      NUMSQ=NUMSQ+1
      SUMSQN=SUMSQN+(DHMF2N-HMF2N)**2
      SUMSQS=SUMSQS+(DHMF2S-HMF2S)**2
      SUMN=SUMN+(DHMF2N-HMF2N)
      SUMS=SUMS+(DHMF2S-HMF2S)

      RETURN
      END
C:::::::::::::::::::::::::: ADJDTMAX :::::::::::::::::::::
C..... This subroutine is used with ExB drift to shorten time step near the end
C..... of the run so that the run ends as close as possible to specified time
C..... Written by P. Richards December 2008
      SUBROUTINE ADJDTMAX(ITAU,     !.. Total run time
     >                      DT,     !.. Current time step
     >                   TSTOP,     !.. MLength of run
     >                   TWIDT,     !.. Twilight time step
     >                   DTMAX)     !.. Maximum time step
      IMPLICIT NONE
      INTEGER ITAU
      DOUBLE PRECISION DT,TSTOP,TWIDT,DTMAX
      !.. Don't allow max time step to be too big for ExB
      IF(DTMAX.GT.3000) THEN
         WRITE(6,*) ' ExB drift: Max time step reset to 3000 secs'
         DTMAX=3000
      ENDIF
      !.. Shorten time step near end of the run to end closer to end L
      IF((REAL(ITAU)+DT)/3600.0.GT.TSTOP-0.5) THEN
        IF(DTMAX.GT.300) DTMAX=300
        IF((REAL(ITAU)+DT)/3600.0.GT.TSTOP-0.1.AND.DTMAX.GT.60) DTMAX=60
        TWIDT=DTMAX   !.. make sure the results are stored in PRINT file
      ENDIF
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C...... This subroutine calculates the limiting H+ flux in the Northern Hemisphere
C...... using the method described by Richards and Torr JGR 1985 page 5261
C...... Updated July 2017
      SUBROUTINE FLXLIMIT(JR,  !.. Index of limiting flux starting altitude. See below
     >                    JQ,  !.. Index of equatorial grid point
     >                    JK,  !.. Index of NmF2 grid point
     >                     Z,  !.. Altitude grid (km)
     >                     N,  !.. Ion densitie array O+, H+, tot+, He+ (cm-3)
     >                    TI,  !.. Ion and electron temperatures (K)
     >                    BM,  !.. Magnetic field strength
     >                    GR,  !.. Gravity array (cm-2s-1)
     >                    ON,  !.. Atomic oxygen density array (cm-3)
     >                    HN,  !.. Atomic hydrogen density array (cm-3)
     >                    TN,  !.. Neutral temperature array (K)
     >                    SL,  !.. Distance along field line (cm)
     >                    JZ,  !.. OUT: Index of ZTRAN
     >                 FYLIM,  !.. OUT: Limiting flux (cm-2 s-1)
     >                HPLOSS,  !.. OUT: H+ total loss above escape altitude (cm-2)
     >                    HS,  !.. OUT: O+ scale height (cm)
     >                 ZTRAN)  !.. OUT: Altitude of transition to escape flux (km)
      INTEGER JR,JQ,JK
      DOUBLE PRECISION RTS(199),ON(401),HN(401),TN(401),BM(401),SL(401)
     > ,Z(401),N(4,401),TI(3,401),GR(2,401),H1,H2,H0,ACON,ZTRAN,HS,
     > FYLIM,HPLOSS
     
      !.. JR is the grid point corresponding to the reference altitude (Zr) 250 to
      !.. 300km below the expected transition altitude. FLIP uses
      !.. Zr = 200.0+200.0*(F107+F107A)/2.0/70.

      !.. Determine the H+ flux by simple method
      H1=-51.6*(TI(2,JR)+TI(3,JR))/GR(2,JR)  !.. O+ scale height
      H2=-51.6*TN(JR)/GR(2,JR)               !.. O scale height
      H0=1.0E5*H1*H2/(H1-H2)
      
      !.. H+ + O and O+ + H reaction rates are from Fox and Sung JGR 2001, page 21,305
      CALL RATS(JR,TI(3,JR),TI(2,JR),TN(JR),RTS)  !.. get reaction rates

      !.. Calculate the transition altitude
      ACON=9E7*TI(2,JR)**2.5/RTS(1)/H0**2
      ZTRAN=Z(JR)-DLOG(ACON/N(1,JR)/ON(JR))*H1*H2/(H1+H2)

      !.. Find the index of the transition altitude index JZ
      JZ=JK   !.. starting at hmF2 altitude and ascending
      DO J=1,JQ
        IF(J.GT.JK.AND.Z(J).LE.ZTRAN) JZ=J
      ENDDO

C       !.. Alternative for finding transition altitude. Use the peak H+ density.
C       It is lower than ZTRAN and overestimates the limiting flux by about 25%
C      JZ=JQ   !.. Begin at equator and go down the field line towards hmF2
C      DO J=JQ-1,JK,-1
C        !.. Check for peak density altitude
C        IF(N(2,J).GT.N(2,J-1).AND.N(2,J).GT.N(2,J+1))
C     >    JZ=J
C      ENDDO

      IF(ABS(Z(JZ)-Z(JR)).GT.500.0) JZ=JR  !.. error capture
      
      HS=-51.6E5*(TI(2,JZ)+TI(3,JZ))/GR(2,JZ)      !.. O+ scale height
      !.. Limiting flux normalized to the first grid point cross section
      !.. For comparison, the FLIP flux is also normalized to the same 
      !.. altitude for output
      FYLIM=RTS(2)*HN(JZ)*N(1,JZ)*HS*BM(1)/BM(JZ)

      !.. Integrate the H+ total loss above ZTRAN. Should be small
      HPLOSS=0.0D0
      DO J=JZ,JQ
        HPLOSS=HPLOSS+RTS(1)*N(2,J)*ON(J)*(SL(J+1)-SL(J))
      ENDDO
      RETURN
      END
C:::::::::::::::::::::::::::::: FLIP_MODS ::::::::::::::::::::::::::::::::::::::::
C......This routine records changes made since May 2017
      SUBROUTINE FLIP_MODS
      IMPLICIT NONE
      WRITE(6,'(//A)') ' ====================================='
      WRITE(6,'(/A,//3X,A/)') ' ..... FLIP model Modifications .....',
     >  'Search for the date in the FLIP model to find the modification'
      WRITE(6,'(A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,9A)') 
     >  ' 2017-05-20: Modified routine ALTWIN to reduce wind speeds'
     > ,' when the input hmF2 is below normal F2 altitudes (<~215km).'
     > ,' This is mainly a problem when there is no F2 peak, which'
     > ,' is known as the G condition on ionograms. It is'
     > ,' fairly common around midday in summer at solar minimum,'
     > ,' when O+ loss rates are high.'
     > ,' FLIP uses a running mean for missing hmF2, but this may not'
     > ,' represent the actual situation. As a result, the hmF2 winds'
     > ,' can become excessive and erratic. The modification reduces'
     > ,' the number and size of excessive winds.'
      WRITE(6,'(A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,9A)') 
     >  ' 2017-07-08: Redid the Richards and Torr JGR 1985 limiting'
     >  ,' flux.  Results are in the end columns of the *Heat*.txt file'
      WRITE(6,'(A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,9A)') 
     >  ' 2017-09-08: Redid solar eclipse routine (ECFUN) to split the'
     >  ,' irradiance into chromosphere responsible for photoinization'
     >  ,' and corona responsible for photoelectrons. Subsequently'
     >  ,' modified to vary 30.4 nm as chromosphere emission'
      WRITE(6,'(A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,9A)') 
     >  ' 2018-02-27: Odd N and N2(v) solns switched off for L < 1.07.'
     >  ,' A warning is provided for ExB drift for L < 1.15'
      WRITE(6,'(A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,9A)') 
     >  ' 2018-03-16: RATCHK routine modified to test reaction rates'
     >  ,' at 3 different temperatures.'
      WRITE(6,'(A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,9A)') 
     >  ' 2018-10-24: Fixed a bug when using a file for EUV.'
     >  ,' Photoionization was not calculated in print phase.'     
      WRITE(6,'(A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,9A)') 
     >  ' 2018-10-30: Updated IRI95 routine for when a requested year'
     >  ,' is greater than the end year in the IG-RZ index file.'    
     >  ,' Now, a message explains how to update the file.'    
      WRITE(6,'(/A,A)') 
     >  ' 2018-11-19: Further updated limiting flux calculation and'
     >  ,' its output to the *HEAT*.txt file.'    
      WRITE(6,'(/A,A)') 
     >  ' 2018-12-03: Made the Kp output real in the *LOG* file.'
     >  ,' Also made it more accurate by using a lookup table'    
      WRITE(6,'(/A,A/2X,A)') 
     >  ' 2020-07-23: Fixed a problem in CALVN2 with N2v factor just'
     >  ,' below the field line grid apex point at low latitudes.'
     >  ,' Also modified vertical grid spacing in VERGRD.'
      WRITE(6,'(/A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,9A)') 
     >  ' 2020-07-23: Created routine TEST_WIND_FILES to test wind/hmF2'
     >  ,' file for winds or hmF2. It is used at low latitudes where'
     >  ,' it is not possible to get winds from hmF2. It allows'
     >  ,' winds to be input from a file'
     >  ,'2020-07-29: North winds are positive in both files now.'
      WRITE(6,'(/A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,A,A,/2X,9A)') 
     >  ' 2020-07-29: Modified routine XION in MINORA.FOR to'
     >  ,' chemically calculate He+ when IHEP=0, and N+ when' 
     >  ,' INPLUS=0 in FLIPRUN.DDD'
      WRITE(6,'(/A,A,A)') 
     >  ' 2020-08-02: Fixed a problem in calculating the index of hmF2'
     >  ,' when there is no peak above 180 km in routine DIAGPR.'
      WRITE(6,'(/A,A,/2X,A)') 
     >  ' 2020-08-02: Modified routine CEMOX to output 6300A and 7755A'
     >  ,' emissions for Bob Kerr. TWI-EMIS was updated to accomodate'
     >  ,' the new FLIP-Print output'
      WRITE(6,'(/1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A)') 
     >  '2020-08-30: Changes for getting output above a near equatorial'
     >  ,'station. Multiple runs are made with different L-values'
     >  ,'to produce an electron density profile. Output is on unit 43.'
     >  ,'It is turned on by specifying magnetic latitude with'
     >  ,'the 8th parameter in FLIPRUN.DDD, which was always zero.'
     >  ,'For backward compatibility, if it is < 0.001, no output.'
      WRITE(6,'(/1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A)') 
     >  '2020-09-10: Improved error messages for inputting winds' 
     >  ,'at high and equatorial latitudes'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A)') 
     >  '2020-09-27: Added HEATFS and TOPTENSV CAL_EHEAT to smooth out'
     >  ,' big jumps in Te data. Also since the Te algorithm needs a Te'
     >  ,' difference to work, a small percentage is added to TOPTEN'
      WRITE(6,'(/1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A)') 
     >  '2020-10-28: Fixed the Ti/Tn bug in Rates.for'
      WRITE(6,'(/1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A)') 
     >  '2020-11-06: Took care of TIMEOF < TIMEON in FLIPRUN.DDD. By'
     >  ,'FLIP adding 24 to TIMEOF for the run.'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A)') 
     >  '2021-02-02: The reduction in plasmaspheric H+ and He+ after'
     >  ,'high Kp was changed to Chappell Rev. Geophys 1972 figure 4.'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A)')
     >  '2021-07-15: Implemented Apex coordinates. FLIP may need'
     >  ,'updating after 2025 but only if north and south geographic'
     >  ,'latitude and longitude are too different from the IGRF'
     >  ,'values. To update, replace the Apex.f90 file in the folder'
     >  ,'FLIP with the latest version from the Web and recompile.'
     >  ,'https://github.com/NCAR/apex_fortran/blob/master/apex.f90'
     >  ,'https://omniweb.gsfc.nasa.gov/vitmo/cgm.html'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A)')
     >  '2022-05-16: Modified the hmF2 proxy wind algorithm to allow'
     >  ,'hmF2 to go down to 195 km based on Kharkiv December 2019'
     >  ,'radar data showing hmF2 at 200 km'
      WRITE(6,'(/1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A)')
     >  '2022-05-29: FLIPRUN.DDD modified to allow different H density'
     >  ,'multipliers for the Northern and Southern hemispheres.'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A)')
     >  '2022-08-06: Fixed a problem with the field line grid for'
     >  ,'low L values when the field line does not reach 250 km.'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A)')
     >  '2022-08-15: Modification to allow Hydrogen densities to be'
     >  ,'read in from the hmF2/wind data files'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A)')
     >  '2022-10-15: RSLPSD and RSTEMC_ExB modified to output Te heat'
     >  ,'equation details for Bill Peterson'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A)')
     >  '2022-11-01: Modification to print file headers in most files'
     >  ,'every 24 hours at midnight LT'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A)')
     >  '2023-01-09: FLIP run & print files modified to write altitude'
     >  ,'profiles of HWM14 (zon. & mer.) winds & FLIP magnetic'
     >  ,' meridional winds. Also added winds at a specified altitude'
     >  ,'to the UT LT SZA ... output line at the FLIP-print stage'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A)')
     >  '2023-02-06: Modified run & print files to process gravity'
     >  ,'wave data from a file'
      WRITE(6,'(1X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A,1X,A,/3X,A)')
     >  '2023-02-16: Changed UT output from F8.2 to F8.3'
      RETURN
      END